Study 1 Analysis

library(psych)
library(ggplot2)
library(tidyverse)
library(effects)
## Library for the multiple comparisons
library(phia)
## Library to compute the effect sizes
library(sjstats)
library(lsr)
library(lavaan)
library(ggthemes)
library(tibble)
library(car)
library(olsrr)
library(pwr)

library(R2jags)
library(rjags)
library(runjags)
library("lattice")
library("superdiag")
library(gtools)

library(bayestestR) 
library("ggmcmc")
library(brms)
library(BayesFactor)
library(MCMCvis)

source('Useful Code V2.R')

library(readr)
attobo <- read_csv("attobo1_clean.csv")

colnames(attobo) <- tolower(colnames(attobo))

attobo$scenario <- as.factor(attobo$scenario)

attobo$attribution <- attobo$scenario

attobo$attribution_effort <- factor(attobo$attribution, levels = c("Effort", "Strategy", "Aptitude"))
attobo$attribution_strat <- factor(attobo$attribution, levels = c( "Strategy","Effort", "Aptitude"))

attobo$expectancy <- (attobo$expectancy - 1)/10
attobo$expectancy <- ifelse(attobo$expectancy == 1, .999,
                     ifelse(attobo$expectancy == 0, .001, attobo$expectancy) )

Descriptive Stats

## Sample Breakdown
group_by(attobo, attribution) %>%
  summarise(n())

## means and sd of DVs
aggregate(cbind(combined_method, bruteforce, seekhelp)~attribution, data=attobo, FUN=mean)

aggregate(cbind(combined_method, bruteforce, seekhelp)~attribution, data=attobo, FUN=sd)

aggregate(cbind(total, convergent_method)~attribution, data=attobo, FUN=mean)

aggregate(cbind(total, convergent_method)~attribution, data=attobo, FUN=sd)

## Locus, Variability, and Uncontrollability
attobo %>%
  group_by(attribution) %>%
  summarise(meanLocusOut = mean(locusout, na.rm = T),
            sdLocusOut = sd(locusout, na.rm = T),
            meanVar = mean(variability, na.rm = T),
            sdVar = sd(variability, na.rm = T),
            meanUncon = mean(uncontrollability, na.rm = T),
            sdUncon = sd(uncontrollability, na.rm = T))

attobo.advice <- attobo %>%
  select(combined_method, bruteforce, seekhelp, convergent_method,
         total, expectancy, total_method)

attobo.attribution <- attobo %>%
  select(aptout, aptvar, aptuncon, 
         effortout, effortvar, effortuncon,
         stratout, stratvar, stratuncon,
         luckout, luckvar, luckuncon,
         taskout, taskvar, taskuncon,
         moodout, moodvar, mooduncon)

attobo.others <- attobo %>%
  select(expectancy, gender, age, education, sms_mean,
         smo_mean, mp_mean, so_mean)
  

ggplot(data = gather(attobo.advice, factor_key = TRUE), aes(x = factor(value))) +
  geom_bar() + 
  facet_wrap(~ key,scales = "free", as.table = TRUE, nrow = 3) + xlab("") + 
  theme_bw()

ggplot(data = gather(attobo.attribution, factor_key = TRUE), aes(x = factor(value))) +
  geom_bar() + 
  facet_wrap(~ key,scales = "free", as.table = TRUE, nrow = 3) + xlab("") + 
  theme_bw()

ggplot(data = gather(attobo.others, factor_key = TRUE), aes(x = factor(value))) +
  geom_bar() + 
  facet_wrap(~ key,scales = "free", as.table = TRUE, nrow = 3) + xlab("") + 
  theme_bw()

Power Analysis

pwr.chisq.test(w = .3, df = 2, sig.level =.05, power = .8)

## 
##      Chi squared power calculation 
## 
##               w = 0.3
##               N = 107.0521
##              df = 2
##       sig.level = 0.05
##           power = 0.8
## 
## NOTE: N is the number of observations

Frequentist Estimate

library(glm2)
library(betareg)

method.freq <- glm2(combined_method ~ attribution_strat, data = attobo, family = poisson(link = "log"))

convergent.freq <- glm2(convergent_method ~ attribution_effort, data = attobo, family = poisson(link = "log"))

bruteforce.freq <- glm2(bruteforce ~ attribution_effort, data = attobo, family = poisson(link = "log"))

expectancy.freq <- betareg(expectancy ~ attribution, data = attobo, link = "logit" )

seekhelp.freq <- glm2(seekhelp ~ attribution, data = attobo, family = poisson(link = "log"))

total.freq <- glm2(total ~ attribution, data = attobo, family = poisson(link = "log"))


summary(glm2(combined_method ~ attribution_strat, data = attobo, family = poisson(link = "log")))

## 
## Call:
## glm2(formula = combined_method ~ attribution_strat, family = poisson(link = "log"), 
##     data = attobo)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.74895  -1.18818   0.09237   0.36315   2.69974  
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                0.42488    0.09806   4.333 1.47e-05 ***
## attribution_stratEffort   -0.77319    0.17449  -4.431 9.38e-06 ***
## attribution_stratAptitude -0.51870    0.16127  -3.216   0.0013 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 209.38  on 202  degrees of freedom
## Residual deviance: 186.37  on 200  degrees of freedom
## AIC: 508.8
## 
## Number of Fisher Scoring iterations: 5

summary(glm2(convergent_method ~ attribution_effort, data = attobo, family = poisson(link = "log")))

## 
## Call:
## glm2(formula = convergent_method ~ attribution_effort, family = poisson(link = "log"), 
##     data = attobo)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4038  -0.8639  -0.5688   0.4310   1.8607  
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                -0.01482    0.12217  -0.121    0.903    
## attribution_effortStrategy -1.80680    0.32528  -5.555 2.78e-08 ***
## attribution_effortAptitude -0.97100    0.23435  -4.143 3.42e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 209.45  on 202  degrees of freedom
## Residual deviance: 161.03  on 200  degrees of freedom
## AIC: 339.06
## 
## Number of Fisher Scoring iterations: 5

summary(glm2(bruteforce ~ attribution_effort, data = attobo, family = poisson(link = "log")))

## 
## Call:
## glm2(formula = bruteforce ~ attribution_effort, family = poisson(link = "log"), 
##     data = attobo)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.8911  -0.6229  -0.6229  -0.3430   2.6321  
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -0.9237     0.1925  -4.800 1.59e-06 ***
## attribution_effortStrategy  -1.9095     0.5358  -3.564 0.000365 ***
## attribution_effortAptitude  -0.7161     0.3376  -2.121 0.033905 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 167.14  on 202  degrees of freedom
## Residual deviance: 147.77  on 200  degrees of freedom
## AIC: 226.05
## 
## Number of Fisher Scoring iterations: 6

summary(betareg(expectancy ~ attribution, data = attobo, link = "logit" ))

## 
## Call:
## betareg(formula = expectancy ~ attribution, data = attobo, link = "logit")
## 
## Standardized weighted residuals 2:
##     Min      1Q  Median      3Q     Max 
## -4.6329 -0.4291 -0.1732  0.1347  3.3119 
## 
## Coefficients (mean model with logit link):
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          0.740427   0.130347   5.680 1.34e-08 ***
## attributionEffort    0.440632   0.181451   2.428   0.0152 *  
## attributionStrategy -0.002054   0.179944  -0.011   0.9909    
## 
## Phi coefficients (precision model with identity link):
##       Estimate Std. Error z value Pr(>|z|)    
## (phi)   2.4961     0.2267   11.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Type of estimator: ML (maximum likelihood)
## Log-likelihood:  64.7 on 4 Df
## Pseudo R-squared: 0.03119
## Number of iterations: 19 (BFGS) + 2 (Fisher scoring)

summary(glm2(total ~ attribution, data = attobo, family = poisson(link = "log")))

## 
## Call:
## glm2(formula = total ~ attribution, family = poisson(link = "log"), 
##     data = attobo)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.22288  -0.75504   0.03143   0.51507   2.34979  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          0.78574    0.08248   9.527   <2e-16 ***
## attributionEffort    0.11872    0.11294   1.051    0.293    
## attributionStrategy -0.11490    0.11967  -0.960    0.337    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 128.65  on 202  degrees of freedom
## Residual deviance: 124.56  on 200  degrees of freedom
## AIC: 650.56
## 
## Number of Fisher Scoring iterations: 4

Anova(glm2(combined_method ~ attribution_strat, data = attobo, family = poisson(link = "log")))

Anova(glm2(convergent_method ~ attribution_effort, data = attobo, family = poisson(link = "log")))

Anova(glm2(bruteforce ~ attribution_effort, data = attobo, family = poisson(link = "log")))

Anova(betareg(expectancy ~ attribution, data = attobo, link = "logit" ))

Anova(glm2(total ~ attribution, data = attobo, family = poisson(link = "log")))

Bayesian Model - method Advice

strategy <- ifelse(attobo$attribution == "Strategy", 1,0)
effort <- ifelse(attobo$attribution == "Effort", 1, 0) 

method <- attobo$combined_method

attobo.datjags <- list(strategy = strategy, effort = effort, method = method, N = nrow(attobo))
attobo.datjags

Model Specification

attobo.model <- function() {
  for(i in 1:N){
    # lambda is mean and variance of the poisson distribution
    method[i] ~ dpois(lambda[i])
    log(lambda[i]) <- b[1] + b[2] * effort[i] + b[3] * strategy[i]
  }
  
  #priors
  b[1] ~ dnorm(0, 0.0025)
  b[2] ~ dnorm(0, 0.0025)
  b[3] ~ dnorm(0, 0.0025)
  
  for (i in 1:N){
    
    error[i] <- lambda[i] - method[i]
  
  }
  
  
}

attobo.params <- c("b", "lambda", "error")

b.inits <- method.freq$coefficients

attobo.inits1 <-  list( "b" = b.inits)
attobo.inits2 <- list( "b" = rep(0, 3) )
attobo.inits3 <- list( "b" = rep(1, 3) )
attobo.inits4 <- list( "b" = rep(-1, 3) )
attobo.inits5 <- list( "b" = rep(2, 3) )
attobo.inits <- list(attobo.inits1, attobo.inits2,attobo.inits3, attobo.inits4,attobo.inits5)

set.seed(1128)
attobo.fit <- jags(data = attobo.datjags, inits = attobo.inits, parameters.to.save = attobo.params, 
                n.chains = 5, n.iter = 100000, n.burnin = 10000, model.file = attobo.model)

attobo.fit.upd <- update(attobo.fit, n.iter =1000)
attobo.fit.upd <- autojags(attobo.fit)

attobo.fit.mcmc <- as.mcmc(attobo.fit)
save(attobo.fit.mcmc, file = "attobo.method.mcmc.RData")

Convergence Diagnostics

load("attobo.method.mcmc.RData")

MCMCtrace(attobo.fit.mcmc[,1:4], ind = TRUE, pdf = FALSE )

autocorr.plot(attobo.fit.mcmc[1][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[2][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[3][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[4][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[5][,1:4], ask = FALSE)

attobo.out <- as.data.frame(as.matrix(attobo.fit.mcmc))
## saving data
save(attobo.out, file = "attobo.method.RData")

load("attobo.method.RData")

pred.method <- attobo.out[, grep("lambda[", colnames(attobo.out), fixed = T)]

## estimated count for each participant
pred.method <- pred.method[, c(mixedsort(names(pred.method)))]

strategy.id <- as.numeric(unlist(which(attobo$attribution == "Strategy")))
effort.id <- as.numeric(unlist(which(attobo$attribution == "Effort")))
aptitude.id <- as.numeric(unlist(which(attobo$attribution == "Aptitude")))

pred.method$strategy.sim <- apply(pred.method[, strategy.id], 1, mean)
pred.method$effort.sim <- apply(pred.method[, effort.id], 1, mean)
pred.method$aptitude.sim <- apply(pred.method[, aptitude.id], 1, mean)

quantile(pred.method$strategy.sim, probs = c(.025, .5, .975))

##     2.5%      50%    97.5% 
## 1.250904 1.526093 1.835312

quantile(pred.method$effort.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.5170229 0.7007571 0.9303260

quantile(pred.method$aptitude.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.6974865 0.9049689 1.1477881

pred.method$diff.1 <- pred.method$strategy.sim - pred.method$aptitude.sim
pred.method$diff.2 <- pred.method$strategy.sim - pred.method$effort.sim
pred.method$diff.3 <- pred.method$effort.sim - pred.method$aptitude.sim

quantile(pred.method$diff.1, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.2514705 0.6191405 0.9899373

quantile(pred.method$diff.2, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.4776900 0.8261376 1.1801479

quantile(pred.method$diff.3, probs = c(.025, .5, .975))

##       2.5%        50%      97.5% 
## -0.5085726 -0.2057398  0.1046953

error.method <- attobo.out[, grep("error[", colnames(attobo.out), fixed = T)]

error <- as.vector(unlist(lapply(error.method, mean)))
mean(error)

## [1] 0.001381898

plot(density(error))

Expected Counts for each attribution

attobo.posterior.coefs <- select(pred.method, `strategy.sim`, `effort.sim`,`aptitude.sim`)
names(attobo.posterior.coefs) <- c("Strategy", "Effort", "Aptitude")

attobo.posterior.coefs.long <- gather(attobo.posterior.coefs)
#head(attobo.posterior.coefs.long)

library("ggridges")
ggplot(data = attobo.posterior.coefs.long, aes(x = value, y = key)) + 
  stat_density_ridges(quantile_lines = TRUE, 
                      quantiles = c(0.025, 0.5, 0.975), 
                      alpha = 0.7)

Bayesian Model - bruteforce Advice

strategy <- ifelse(attobo$attribution == "Strategy", 1,0)
effort <- ifelse(attobo$attribution == "Effort", 1, 0) 

bruteforce <- attobo$bruteforce

attobo.datjags <- list(strategy = strategy, effort = effort, bruteforce = bruteforce, N = nrow(attobo))
attobo.datjags

Model Specification

attobo.model <- function() {
  for(i in 1:N){
    # lambda is mean and variance of the poisson distribution
    bruteforce[i] ~ dpois(lambda[i])
    log(lambda[i]) <- b[1] + b[2] * effort[i] + b[3] * strategy[i]
  }
  
  #priors
  b[1] ~ dnorm(0, 0.0025)
  b[2] ~ dnorm(0, 0.0025)
  b[3] ~ dnorm(0, 0.0025)
  
  for (i in 1:N){
    
    error[i] <- lambda[i] - bruteforce[i]
  
  }
  
  
}

attobo.params <- c("b", "lambda", "error")

b.inits <- bruteforce.freq$coefficients

attobo.inits1 <-  list( "b" = b.inits)
attobo.inits2 <- list( "b" = rep(0, 3) )
attobo.inits3 <- list( "b" = rep(1, 3) )
attobo.inits4 <- list( "b" = rep(-1, 3) )
attobo.inits5 <- list( "b" = rep(2, 3) )
attobo.inits <- list(attobo.inits1, attobo.inits2,attobo.inits3, attobo.inits4,attobo.inits5)

set.seed(1128)
attobo.fit <- jags(data = attobo.datjags, inits = attobo.inits, parameters.to.save = attobo.params, 
                n.chains = 5, n.iter = 100000, n.burnin = 10000, model.file = attobo.model)

attobo.fit.upd <- update(attobo.fit, n.iter =1000)
attobo.fit.upd <- autojags(attobo.fit)

attobo.fit.mcmc <- as.mcmc(attobo.fit)
save(attobo.fit.mcmc, file = "attobo.bruteforce.mcmc.RData")

Convergence Diagnostics

load("attobo.bruteforce.mcmc.RData")

MCMCtrace(attobo.fit.mcmc[,1:4], ind = TRUE, pdf = FALSE )

autocorr.plot(attobo.fit.mcmc[1][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[2][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[3][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[4][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[5][,1:4], ask = FALSE)

attobo.out <- as.data.frame(as.matrix(attobo.fit.mcmc))

## saving data
save(attobo.out, file = "attobo.bruteforce.RData")

pred.bruteforce <- attobo.out[, grep("lambda[", colnames(attobo.out), fixed = T)]

## estimated count for each participant
pred.bruteforce <- pred.bruteforce[, c(mixedsort(names(pred.bruteforce)))]

strategy.id <- as.numeric(unlist(which(attobo$attribution == "Strategy")))
effort.id <- as.numeric(unlist(which(attobo$attribution == "Effort")))
aptitude.id <- as.numeric(unlist(which(attobo$attribution == "Aptitude")))

pred.bruteforce$strategy.sim <- apply(pred.bruteforce[, strategy.id], 1, mean)
pred.bruteforce$effort.sim <- apply(pred.bruteforce[, effort.id], 1, mean)
pred.bruteforce$aptitude.sim <- apply(pred.bruteforce[, aptitude.id], 1, mean)

quantile(pred.bruteforce$strategy.sim, probs = c(.025, .5, .975))

##     2.5%      50%    97.5% 
## 1.250904 1.526093 1.835312

quantile(pred.bruteforce$effort.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.5170229 0.7007571 0.9303260

quantile(pred.bruteforce$aptitude.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.6974865 0.9049689 1.1477881

error.bruteforce <- attobo.out[, grep("error[", colnames(attobo.out), fixed = T)]

error <- as.vector(unlist(lapply(error.bruteforce, mean)))
mean(error)

## [1] 0.001381898

plot(density(error))

Expected Counts for each attribution

attobo.posterior.coefs <- select(pred.bruteforce, `strategy.sim`, `effort.sim`,`aptitude.sim`)
names(attobo.posterior.coefs) <- c("Strategy", "Effort", "Aptitude")

attobo.posterior.coefs.long <- gather(attobo.posterior.coefs)
#head(attobo.posterior.coefs.long)

library("ggridges")
ggplot(data = attobo.posterior.coefs.long, aes(x = value, y = key)) + 
  stat_density_ridges(quantile_lines = TRUE, 
                      quantiles = c(0.025, 0.5, 0.975), 
                      alpha = 0.7)

Bayesian Model - mandatory Advice

strategy <- ifelse(attobo$attribution == "Strategy", 1,0)
effort <- ifelse(attobo$attribution == "Effort", 1, 0) 

convergent <- attobo$convergent_method

attobo.datjags <- list(strategy = strategy, effort = effort, convergent = convergent, N = nrow(attobo))
attobo.datjags

Model Specification

attobo.model <- function() {
  for(i in 1:N){
    # lambda is mean and variance of the poisson distribution
    convergent[i] ~ dpois(lambda[i])
    log(lambda[i]) <- b[1] + b[2] * effort[i] + b[3] * strategy[i]
  }
  
  #priors
  b[1] ~ dnorm(0, 0.0025)
  b[2] ~ dnorm(0, 0.0025)
  b[3] ~ dnorm(0, 0.0025)
  
  for (i in 1:N){
    
    error[i] <- lambda[i] - convergent[i]
  
  }
  
  
}

attobo.params <- c("b", "lambda", "error")

b.inits <- convergent.freq$coefficients

attobo.inits1 <-  list( "b" = b.inits)
attobo.inits2 <- list( "b" = rep(0, 3) )
attobo.inits3 <- list( "b" = rep(1, 3) )
attobo.inits4 <- list( "b" = rep(-1, 3) )
attobo.inits5 <- list( "b" = rep(2, 3) )
attobo.inits <- list(attobo.inits1, attobo.inits2,attobo.inits3, attobo.inits4,attobo.inits5)

set.seed(1128)
attobo.fit <- jags(data = attobo.datjags, inits = attobo.inits, parameters.to.save = attobo.params, 
                n.chains = 5, n.iter = 100000, n.burnin = 10000, model.file = attobo.model)

attobo.fit.upd <- update(attobo.fit, n.iter =1000)
attobo.fit.upd <- autojags(attobo.fit)

attobo.fit.mcmc <- as.mcmc(attobo.fit)
save(attobo.fit.mcmc, file = "attobo.convergent.mcmc.RData")

Convergence Diagnostics

load("attobo.convergent.mcmc.RData")
MCMCtrace(attobo.fit.mcmc[,1:4], ind = TRUE, pdf = FALSE )

autocorr.plot(attobo.fit.mcmc[1][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[2][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[3][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[4][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[5][,1:4], ask = FALSE)

attobo.out <- as.data.frame(as.matrix(attobo.fit.mcmc))

## saving data
save(attobo.out, file = "attobo.convergent.RData")

load("attobo.convergent.RData")
pred.convergent <- attobo.out[, grep("lambda[", colnames(attobo.out), fixed = T)]

## estimated count for each participant
pred.convergent <- pred.convergent[, c(mixedsort(names(pred.convergent)))]

strategy.id <- as.numeric(unlist(which(attobo$attribution == "Strategy")))
effort.id <- as.numeric(unlist(which(attobo$attribution == "Effort")))
aptitude.id <- as.numeric(unlist(which(attobo$attribution == "Aptitude")))

pred.convergent$strategy.sim <- apply(pred.convergent[, strategy.id], 1, mean)
pred.convergent$effort.sim <- apply(pred.convergent[, effort.id], 1, mean)
pred.convergent$aptitude.sim <- apply(pred.convergent[, aptitude.id], 1, mean)

quantile(pred.convergent$strategy.sim, probs = c(.025, .5, .975))

##       2.5%        50%      97.5% 
## 0.08305286 0.15728183 0.26887682

quantile(pred.convergent$effort.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.7675159 0.9802979 1.2365364

quantile(pred.convergent$aptitude.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.2390818 0.3637082 0.5336015

error.convergent <- attobo.out[, grep("error[", colnames(attobo.out), fixed = T)]

error <- as.vector(unlist(lapply(error.convergent, mean)))
mean(error)

## [1] 0.000883394

plot(density(error))

Expected Counts for each attribution

attobo.posterior.coefs <- select(pred.convergent, `strategy.sim`, `effort.sim`,`aptitude.sim`)
names(attobo.posterior.coefs) <- c("Strategy", "Effort", "Aptitude")

attobo.posterior.coefs.long <- gather(attobo.posterior.coefs)
#head(attobo.posterior.coefs.long)

library("ggridges")
ggplot(data = attobo.posterior.coefs.long, aes(x = value, y = key)) + 
  stat_density_ridges(quantile_lines = TRUE, 
                      quantiles = c(0.025, 0.5, 0.975), 
                      alpha = 0.7)

Bayesian Model - seekhelp Advice

strategy <- ifelse(attobo$attribution == "Strategy", 1,0)
effort <- ifelse(attobo$attribution == "Effort", 1, 0) 

seekhelp <- attobo$seekhelp

attobo.datjags <- list(strategy = strategy, effort = effort, seekhelp = seekhelp, N = nrow(attobo))
attobo.datjags

Model Specification

attobo.model <- function() {
  for(i in 1:N){
    # lambda is mean and variance of the poisson distribution
    seekhelp[i] ~ dpois(lambda[i])
    log(lambda[i]) <- b[1] + b[2] * effort[i] + b[3] * strategy[i]
  }
  
  #priors
  b[1] ~ dnorm(0, 0.0025)
  b[2] ~ dnorm(0, 0.0025)
  b[3] ~ dnorm(0, 0.0025)
  
  for (i in 1:N){
    
    error[i] <- lambda[i] - seekhelp[i]
  
  }
  
  
}

attobo.params <- c("b", "lambda", "error")

b.inits <- seekhelp.freq$coefficients

attobo.inits1 <-  list( "b" = b.inits)
attobo.inits2 <- list( "b" = rep(0, 3) )
attobo.inits3 <- list( "b" = rep(1, 3) )
attobo.inits4 <- list( "b" = rep(-1, 3) )
attobo.inits5 <- list( "b" = rep(2, 3) )
attobo.inits <- list(attobo.inits1, attobo.inits2,attobo.inits3, attobo.inits4,attobo.inits5)

set.seed(1128)
attobo.fit <- jags(data = attobo.datjags, inits = attobo.inits, parameters.to.save = attobo.params, 
                n.chains = 5, n.iter = 100000, n.burnin = 10000, model.file = attobo.model)

attobo.fit.upd <- update(attobo.fit, n.iter =1000)
attobo.fit.upd <- autojags(attobo.fit)

attobo.fit.mcmc <- as.mcmc(attobo.fit)
save(attobo.fit.mcmc, file = "attobo.seekhelp.mcmc.RData")

Convergence Diagnostics

load("attobo.seekhelp.mcmc.RData")

MCMCtrace(attobo.fit.mcmc[,1:4], ind = TRUE, pdf = FALSE )

autocorr.plot(attobo.fit.mcmc[1][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[2][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[3][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[4][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[5][,1:4], ask = FALSE)

attobo.out <- as.data.frame(as.matrix(attobo.fit.mcmc))

## saving data
save(attobo.out, file = "attobo.seekhelp.RData")

load("attobo.seekhelp.RData")
pred.seekhelp <- attobo.out[, grep("lambda[", colnames(attobo.out), fixed = T)]

## estimated count for each participant
pred.seekhelp <- pred.seekhelp[, c(mixedsort(names(pred.seekhelp)))]

strategy.id <- as.numeric(unlist(which(attobo$attribution == "Strategy")))
effort.id <- as.numeric(unlist(which(attobo$attribution == "Effort")))
aptitude.id <- as.numeric(unlist(which(attobo$attribution == "Aptitude")))

pred.seekhelp$strategy.sim <- apply(pred.seekhelp[, strategy.id], 1, mean)
pred.seekhelp$effort.sim <- apply(pred.seekhelp[, effort.id], 1, mean)
pred.seekhelp$aptitude.sim <- apply(pred.seekhelp[, aptitude.id], 1, mean)

mean(pred.seekhelp$strategy.sim)

## [1] 0.2078543

quantile(pred.seekhelp$strategy.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.1151397 0.2016900 0.3276598

mean(pred.seekhelp$effort.sim)

## [1] 0.3844671

quantile(pred.seekhelp$effort.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.2542277 0.3784910 0.5406089

mean(pred.seekhelp$aptitude.sim)

## [1] 0.7167676

quantile(pred.seekhelp$aptitude.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.5287659 0.7127448 0.9281304

error.seekhelp <- attobo.out[, grep("error[", colnames(attobo.out), fixed = T)]

error <- as.vector(unlist(lapply(error.seekhelp, mean)))

mean(error)

## [1] 0.001484177

plot(density(error))

Expected Counts for each attribution

attobo.posterior.coefs <- select(pred.seekhelp, `strategy.sim`, `effort.sim`,`aptitude.sim`)
names(attobo.posterior.coefs) <- c("Strategy", "Effort", "Aptitude")

attobo.posterior.coefs.long <- gather(attobo.posterior.coefs)
#head(attobo.posterior.coefs.long)

library("ggridges")
ggplot(data = attobo.posterior.coefs.long, aes(x = value, y = key)) + 
  stat_density_ridges(quantile_lines = TRUE, 
                      quantiles = c(0.025, 0.5, 0.975), 
                      alpha = 0.7)

Expectancy

Bayesian Model - expectancy

strategy <- ifelse(attobo$attribution == "Strategy", 1,0)
effort <- ifelse(attobo$attribution == "Effort", 1, 0) 

attobo.datjags <- list(strategy = strategy, effort = effort, N = nrow(attobo), expectancy = attobo$expectancy)
attobo.datjags

Model Specification

attobo.model <- function() {
  
  for(i in 1:N){
  
    expectancy[i] ~ dbeta(alpha[i], beta[i])
    alpha[i] = phi * mu[i]
    beta[i] = phi * (1 - mu[i])
    logit(mu[i]) <- b[1] + b[2] * effort[i] + b[3] * strategy[i]
  }
  
  #priors
  
 for(i in 1:3){
   b[i] ~ dnorm(0, 0.0025) 
 }
  
    phiinv ~ dgamma(0.01,0.01)
   phi <- 1/phiinv 

  
  #calculations
  for (i in 1:N){
    
    error[i] <- mu[i] - expectancy[i]
    
  }
  
}

attobo.params <- c("b", "mu", "error")

b.inits <- expectancy.freq$coefficients$mean

attobo.inits1 <-  list( "b" = b.inits, "phiinv" = 1)
attobo.inits2 <- list( "b" = rep(0, 3) , "phiinv" =2 )
attobo.inits3 <- list( "b" = rep(1, 3), "phiinv" = 1 )
attobo.inits4 <- list( "b" = rep(-1, 3) , "phiinv" = 2 )
attobo.inits5 <- list( "b" = rep(2, 3), "phiinv" = 1 )
attobo.inits <- list(attobo.inits1, attobo.inits2,attobo.inits3, attobo.inits4,attobo.inits5)

set.seed(1128)
attobo.fit <- jags(data = attobo.datjags, inits = attobo.inits, parameters.to.save = attobo.params, 
                n.chains = 5, n.iter = 100000, n.burnin = 10000, model.file = attobo.model)

attobo.fit.upd <- update(attobo.fit, n.iter =1000)
attobo.fit.upd <- autojags(attobo.fit)

attobo.fit.mcmc <- as.mcmc(attobo.fit)
save(attobo.fit.mcmc, file = "attobo.expectancy.mcmc.RData")

Convergence Diagnostics

load("attobo.expectancy.mcmc.RData")
MCMCtrace(attobo.fit.mcmc[, 1:4] , ind = TRUE, pdf = FALSE )

autocorr.plot(attobo.fit.mcmc[1][, 1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[2][, 1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[3][, 1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[4][, 1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[5][, 1:4], ask = FALSE)

attobo.out <- as.data.frame(as.matrix(attobo.fit.mcmc))

## saving data
save(attobo.out, file = "attobo.expectancy.RData")

load("attobo.expectancy.RData")
pred.expectancy <- attobo.out[, grep("mu[", colnames(attobo.out), fixed = T)]

## estimated count for each participant
pred.expectancy <- pred.expectancy[, c(mixedsort(names(pred.expectancy)))]

strategy.id <- as.numeric(unlist(which(attobo$attribution == "Strategy")))
effort.id <- as.numeric(unlist(which(attobo$attribution == "Effort")))
aptitude.id <- as.numeric(unlist(which(attobo$attribution == "Aptitude")))

pred.expectancy$strategy.sim <- apply(pred.expectancy[, strategy.id], 1, mean)
pred.expectancy$effort.sim <- apply(pred.expectancy[, effort.id], 1, mean)
pred.expectancy$aptitude.sim <- apply(pred.expectancy[, aptitude.id], 1, mean)

quantile(pred.expectancy$strategy.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.6171928 0.6773310 0.7289890

quantile(pred.expectancy$effort.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.7151146 0.7650822 0.8095703

quantile(pred.expectancy$aptitude.sim, probs = c(.025, .5, .975))

##      2.5%       50%     97.5% 
## 0.6205237 0.6769050 0.7312879

mean(pred.expectancy$effort.sim)

## [1] 0.7642647

mean(pred.expectancy$strategy.sim)

## [1] 0.6762033

mean(pred.expectancy$aptitude.sim)

## [1] 0.6765108

error.expectancy <- attobo.out[, grep("error[", colnames(attobo.out), fixed = T)]

error <- as.vector(unlist(lapply(error.expectancy, mean)))
mean(error)

## [1] -0.01727069

plot(density(error))

Expected Counts for each attribution

attobo.posterior.coefs <- select(pred.expectancy, `strategy.sim`, `effort.sim`,`aptitude.sim`)
names(attobo.posterior.coefs) <- c("Strategy", "Effort", "Aptitude")

attobo.posterior.coefs.long <- gather(attobo.posterior.coefs)
#head(attobo.posterior.coefs.long)

library("ggridges")
ggplot(data = attobo.posterior.coefs.long, aes(x = value, y = key)) + 
  stat_density_ridges(quantile_lines = TRUE, 
                      quantiles = c(0.025, 0.5, 0.975), 
                      alpha = 0.7)

Bayesian Model - Total Advice

strategy <- ifelse(attobo$attribution == "Strategy", 1,0)
effort <- ifelse(attobo$attribution == "Effort", 1, 0) 

total <- attobo$total

attobo.datjags <- list(strategy = strategy, effort = effort, total = total, N = nrow(attobo))
attobo.datjags

Model Specification

attobo.model <- function() {
  for(i in 1:N){
    # lambda is mean and variance of the poisson distribution
    total[i] ~ dpois(lambda[i])
    log(lambda[i]) <- b[1] + b[2] * effort[i] + b[3] * strategy[i]
  }
  
  #priors
  b[1] ~ dnorm(0, 0.0025)
  b[2] ~ dnorm(0, 0.0025)
  b[3] ~ dnorm(0, 0.0025)
  
  for (i in 1:N){
    
    error[i] <- lambda[i] - total[i]
  
  }
  
  
}

attobo.params <- c("b", "lambda", "error")

b.inits <- total.freq$coefficients

attobo.inits1 <-  list( "b" = b.inits)
attobo.inits2 <- list( "b" = rep(0, 3) )
attobo.inits3 <- list( "b" = rep(1, 3) )
attobo.inits4 <- list( "b" = rep(-1, 3) )
attobo.inits5 <- list( "b" = rep(2, 3) )
attobo.inits <- list(attobo.inits1, attobo.inits2,attobo.inits3, attobo.inits4,attobo.inits5)

set.seed(1128)
attobo.fit <- jags(data = attobo.datjags, inits = attobo.inits, parameters.to.save = attobo.params, 
                n.chains = 5, n.iter = 100000, n.burnin = 10000, model.file = attobo.model)

attobo.fit.upd <- update(attobo.fit, n.iter =1000)
attobo.fit.upd <- autojags(attobo.fit)

attobo.fit.mcmc <- as.mcmc(attobo.fit)
save(attobo.fit.mcmc, file = "attobo.total.mcmc.RData")

Convergence Diagnostics

load("attobo.total.mcmc.RData")
MCMCtrace(attobo.fit.mcmc[,1:4], ind = TRUE, pdf = FALSE )

autocorr.plot(attobo.fit.mcmc[1][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[2][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[3][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[4][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[5][,1:4], ask = FALSE)

attobo.out <- as.data.frame(as.matrix(attobo.fit.mcmc))

## saving data
save(attobo.out, file = "attobo.total.RData")

load("attobo.total.RData")
pred.total <- attobo.out[, grep("lambda[", colnames(attobo.out), fixed = T)]

## estimated count for each participant
pred.total <- pred.total[, c(mixedsort(names(pred.total)))]

strategy.id <- as.numeric(unlist(which(attobo$attribution == "Strategy")))
effort.id <- as.numeric(unlist(which(attobo$attribution == "Effort")))
aptitude.id <- as.numeric(unlist(which(attobo$attribution == "Aptitude")))

pred.total$strategy.sim <- apply(pred.total[, strategy.id], 1, mean)
pred.total$effort.sim <- apply(pred.total[, effort.id], 1, mean)
pred.total$aptitude.sim <- apply(pred.total[, aptitude.id], 1, mean)

quantile(pred.total$strategy.sim, probs = c(.025, .5, .975))

##     2.5%      50%    97.5% 
## 1.642788 1.950798 2.307109

quantile(pred.total$effort.sim, probs = c(.025, .5, .975))

##     2.5%      50%    97.5% 
## 2.109232 2.462346 2.855788

quantile(pred.total$aptitude.sim, probs = c(.025, .5, .975))

##     2.5%      50%    97.5% 
## 1.857219 2.189896 2.550326

pred.total$diff.1 <- pred.total$strategy.sim - pred.total$aptitude.sim
pred.total$diff.2 <- pred.total$strategy.sim - pred.total$effort.sim
pred.total$diff.3 <- pred.total$effort.sim - pred.total$aptitude.sim

quantile(pred.total$diff.1, probs = c(.025, .5, .975))

##       2.5%        50%      97.5% 
## -0.7115193 -0.2327644  0.2565646

quantile(pred.total$diff.2, probs = c(.025, .5, .975))

##         2.5%          50%        97.5% 
## -1.008393971 -0.510975172 -0.002368577

quantile(pred.total$diff.3, probs = c(.025, .5, .975))

##       2.5%        50%      97.5% 
## -0.2380388  0.2741762  0.8017476

error.total <- attobo.out[, grep("error[", colnames(attobo.out), fixed = T)]

error <- as.vector(unlist(lapply(error.total, mean)))
mean(error)

## [1] 0.001559535

plot(density(error))

Expected Counts for each attribution

attobo.posterior.coefs <- select(pred.total, `strategy.sim`, `effort.sim`,`aptitude.sim`)
names(attobo.posterior.coefs) <- c("Strategy", "Effort", "Aptitude")

attobo.posterior.coefs.long <- gather(attobo.posterior.coefs)
#head(attobo.posterior.coefs.long)

library("ggridges")
ggplot(data = attobo.posterior.coefs.long, aes(x = value, y = key)) + 
  stat_density_ridges(quantile_lines = TRUE, 
                      quantiles = c(0.025, 0.5, 0.975), 
                      alpha = 0.7)